!
!     DARK     Dark frame subtract utility.
!     Copyright (C) 1997-9  Filip Hroch, Masaryk University, Brno, CZ
!

Subroutine Dark(sci, drk, out, ndark)

  implicit none

character(len=*) :: sci, drk, out
integer :: ndark

!      Function Dark(sci1, drk1, out1, nlen)
!      Integer :: Dark, nlen
!      Character :: sci1(nlen), drk1(nlen), out1(nlen)
      
! Parameters

! No. of image dimensions
      Integer, parameter :: DIM = 2

! The darkbat assumes that output image has in intensity
! within range minvalue < x < maxvalue.
! The minvalue always is 0 for measured intensities (?)
! The maxvalue is determined as (2**bitpix - 1) or by SATURATE
! keyword in image headres (in specified order). 
! The specific characteristic of dark images may require some modification
! for real-type intensity of image.
      Real, parameter :: minvalue = 0.0
      Real            :: maxvalue = 65535.0

! FITS variables for darkimage
      Integer :: dstat,dBlockSize,dNaxis,dNaxes(DIM),dPCount,dGCount,dBitPix
      Logical :: dSimple,dExtend
      Real, allocatable, save :: DARKDATA(:,:)

! FITS variables for scientific exposition
      Integer :: stat,sBlockSize,sNaxis,sNaxes(DIM),sPCount,sGCount,BitPix
      Logical :: sSimple,sExtend
      Real, allocatable, save :: SCIDATA(:,:)
      
! FITS variables for output image
      Integer :: ostat
      Real, allocatable, save :: OUTDATA(:,:)

      Logical :: exist, anyf
      Integer :: i
      Real :: x
!      Character(len=256) :: sci, drk, out
      Character(len=256) ::  buf, scalecom, zerocom

!----- 

! C -> Fortran characters conversion
!      do i = 1,nlen
!        sci(i:i) = sci1(i)
!        drk(i:i) = drk1(i) 
!        out(i:i) = out1(i) 
!      enddo

!      write(*,*) '(',trim(sci),') - (',trim(drk),') = (',trim(out),')' 

      dstat = 0
      stat = 0
      ostat = 0

! darkframe initialisation - first image
      if( .not. allocated(DARKDATA) ) then
        Call ftopen(20,drk,0,dBlockSize,dstat)
        Call ftghpr(20,DIM,dSimple,dBitPix,dNaxis,dNaxes,dPCount,dGCount,dExtend,dstat)
        if( dstat /= 0 ) goto 666
        Allocate( DARKDATA(dNaxes(1),dNaxes(2)),stat=dstat)
        Call ftg2de(20,1,minvalue,dNaxes(1),dNaxes(1),dNaxes(2),DARKDATA,anyf,dstat)
        Call ftclos(20,dstat)  
        if( dstat /= 0 ) goto 666
      endif

! scientific exposure initialisation
      Call ftopen(25,sci,0,sBlockSize,stat)
      Call ftghpr(25,DIM,sSimple,BitPix,sNaxis,sNaxes,sPCount,sGCount,sExtend,stat)
      if( stat /= 0 ) goto 666
      if( .not. allocated(SCIDATA) ) &
           Allocate( SCIDATA(sNaxes(1),sNaxes(2)),stat=stat)
      Call ftg2de(25,1,minvalue,sNaxes(1),sNaxes(1),sNaxes(2),SCIDATA,anyf,stat)
      if( stat /= 0 ) goto 666

! a range of intensities of the image
     maxvalue = 2.0**BitPix - 1.0
     call ftgkye(25,'SATURATE',x,buf,stat)
     if( stat == 0 ) maxvalue = x

! output initialisation - parameters by scientific image
      Inquire(file=out,exist=exist,err=666)
      if( exist )then
        open(26,file=out,err=666)
        close(26,status='delete',err=666)
      endif

      if( .not. allocated(OUTDATA) ) &
           Allocate( OUTDATA(sNaxes(1),sNaxes(2)),stat=ostat)
      Call ftinit(26,out,sBlockSize,ostat)
! header keywords copy - by scientific image
      Call ftcopy(25,26,1,ostat)
      Call ftphis(26,'Darkframe subtracted.',ostat)
      if( ostat /= 0 ) goto 666

! dark subtraction
      OUTDATA = SCIDATA - DARKDATA

! range checking
      where( OUTDATA < minvalue )
        OUTDATA = minvalue
      endwhere
      where( OUTDATA > maxvalue )
        OUTDATA = maxvalue
      endwhere

! corrected image output
      call ftp2de(26,1,sNaxes(1),sNaxes(1),sNaxes(2),OUTDATA,ostat)

666   continue

! close files
      if( ostat == 0 ) then     ! delete output (if any error)
        Call ftclos(26,ostat)
      else
        Call ftdelt(26,ostat) 
      endif
      Call ftclos(25,stat)  

      deallocate(OUTDATA,SCIDATA,DARKDATA,stat=ostat)

      ndark = ostat + 1000*stat + 1000000*dstat
      if( (dstat + stat + ostat) /= 0 ) then
        open(3,file=trim(sci)//'.err')
        write(3,*) 'The following error(s) has been occured (',ndark,'):'
        do
          Call ftgmsg(buf)
          if( buf == ' ' ) exit
          write(3,*) trim(buf)
        enddo
        close(3)
      endif


      End
